\documentclass{ctexart} 
\usepackage{amsmath} 
\usepackage{amssymb}
\usepackage{enumerate} 
\usepackage{graphicx} 
\usepackage{subfig}
\usepackage{indentfirst}%首行缩进需要
\usepackage{float}%使图片紧跟文字后面
\CTEXsetup[format = {\Large\bfseries}]{section} \setlength{\parindent}{0pt}
\title{第三次项目实验报告} 
\author{付临\\3200104960\\信息与计算科学} 
\begin{document} 
\maketitle
\setlength{\parindent}{2em}
\section{编译说明}
本次作业采用Makefile对编译进行管理。在终端输入make即可完成编译得到report.pdf文件
，输入make clear可以删除编译过程得到的相应数据以及过程文件。

\section{程序设计说明}
头文件spline.h中设计了函数类Function和插值类Interpolation，并以此生成作业需要是ppForm插值和B样条插值，
同时每种插值方式又根据不同的边界条件分为三种，其中way=1时为complete cubic spline，way=2时为
cubic spline with specified second derivatives at its end points，way=3时为natural
cubic spline。下面将说明每种基类和其派生类的设计思路。
\subsection{Function}
Class Function中设计了operator、diff以及diff\_2，分别返回函数的值、一阶导数和二阶导数。根据两种插值
的需要，派生出class Polynomial和class B\_spline分别返回函数多项式和B样条基函数。
\subsubsection{class Polynomial}
通过ppForm插值将得到分段多项式，因此需要存储每个区间对应的多项式函数，因此设计出class Polynomial，
输入系数$\{c_i\}_i^n$与$x_0$，存储相对应的函数$f(x) = c_{0} + c_{1}(x-x_0) + \cdots +c_{n}(x-x_0)^{n}$。
再通过operoter可以返回分段函数的值。
\subsubsection{class B\_spline}
通过B样条插值需要B样条基函数，输入阶数order和插值点\{t\}即可得到相应阶数和对应区间的函数。根据
Example 3.31可以直接得出order=1时的函数表达式：
\begin{equation*}
    B_i^1(x)=\left\{\begin{aligned}
        &\frac{x-t_{i-1}}{t_{i}-t_{i-1}} &x\in(t_{i-1},t_i]\\
        &\frac{t_{i+1}-x}{t_{i+1}-t_i} &x\in(t_i,t_{i+1}]\\
        &0 &otherwise\\
    \end{aligned}\right.
\end{equation*}
再利用Definition 3.23可以递归得到相应阶数的B样条：
\begin{equation*}
    B_i^{n+1}(x)=\frac{x-t_{i-1}}{t_{i+n}-t_{i-1}}B_i^n(x)
    +\frac{t_{i+n+1}-x}{t_{i+n+1}-t_i}B_{i+1}^n(x)
\end{equation*}
\subsection{Interpolation}
Class Interpolation中设计了插值求解的solve和返回插值函数值的operator，并派生出class
ppForm\_interpolation和B\_spline\_Interpolation来完成插值过程。
\subsubsection{class ppForm\_interpolation}
本次项目作业设计的ppForm插值支持线性插值（order=1）和三次样条插值（order=3），其中三次样条插值
又根据边界条件不同分为三种。输入函数f，插值点、\{t\}、边界条件way、阶数order。

当order=1时，即线性插值。设插值函数为$P(x)$，记$p_i(x)=P(x)|_{[t_i,t_{i+1}]}$，根据
公式：
\begin{equation*}
    p_i(x)=\frac{f(t_{i+1})-f(t_i)}{t_{i+1}-t_i}x+f(t_{i+1})-t_{i}\frac{f(t_{i+1})-f(t_i)}{t_{i+1}-t_i}
\end{equation*}
\noindent 通过solve函数可以直接求解。

当order=3时，即三次样条插值。Solve函数利用Eigen分解LU矩阵可以求解$Ax=b$问题的线性方程组，
从而可以求出插值函数。设插值点为$\{t_i\}_1^n$,根据不同的边界条件，可以得到不同的A和b：\\
$\bullet$way=1
\begin{equation*}
    A=\left(\begin{array}{ccccc}
        1 & & & &\\
        \lambda_2 & 2 & \mu_2 & &\\
        & \ddots & \ddots & \ddots &\\
        & & \lambda_{n-1} & 2 & \mu_{n-1}\\
        & & & & 1\\
    \end{array}\right)
    b=\left(\begin{array}{c}
        f'(t_{1})\\
        3\lambda_2f[t_1,t_2]+3\mu_2f[t_2,t_3]\\
        \vdots \\
        3\lambda_{n-1}f[t_{n-2},t_{n-1}]+3\mu_{n-1}f[t_{n-1},t_{n}]\\
        f'(t_n)\\
    \end{array}\right)
\end{equation*}\\
$\bullet$way=2
\begin{equation*}
    A=\left(\begin{array}{ccccc}
        4 & 2 & & & \\
        \lambda_2 & 2 & \mu_2 & &\\
        & \ddots & \ddots & \ddots &\\
        & & \lambda_{n-1} & 2 & \mu_{n-1}\\
        & & & 2 & 4\\
    \end{array}\right)
    b=\left(\begin{array}{c}
        6f[t_1,t_2]-f''(t_1)(t_2-t_1)\\
        3\lambda_2f[t_1,t_2]+3\mu_2f[t_2,t_3]\\ \vdots \\
        3\lambda_{n-1}f[t_{n-2},t_{n-1}]+3\mu_{n-1}f[t_{n-1},t_{n}]\\
        6f[t_{n-1},t_{n}]+f''(t_n)(t_n-t_{n-1})\\
    \end{array}\right)
\end{equation*}\\
$\bullet$way=3
\begin{equation*}
    A=\left(\begin{array}{ccccc} 
        4 & 2 & & & \\ 
        \lambda_2 & 2 & \mu_2 & &\\
        & \ddots & \ddots & \ddots &\\ 
        & & \lambda_{n-1} & 2 & \mu_{n-1}\\ 
        & & & 2 &4\\ 
    \end{array}\right)
    b=\left(\begin{array}{c}
        6f[t_1,t_2]\\
        3\lambda_2f[t_1,t_2]+3\mu_2f[t_2,t_3]\\ \vdots \\
        3\lambda_{n-1}f[t_{n-2},t_{n-1}]+3\mu_{n-1}f[t_{n-1},t_{n}]\\
        6f[t_{n-1},t_{n}]\\
    \end{array}\right)
\end{equation*}
\noindent 因此通过求解$Ax=b$可以确定唯一的插值多项式，再根据operator可以返回插值函数在任一点
的值。
\subsubsection{class B\_spline\_interpolation}
对于B样条插值，需要加入额外点，故采用等距添加，在插值点两端分别再添加order个额外点。本次作业主要
针对线性插值和三次样条插值。

当order=1时，即线性插值。需要额外插入$t_0$和$t_n$，$t_0=t_1-(t_2-t_1),t_{n+1}=t_n+(t_2-t_1)$。
求解：
\begin{align*}
    Ax&=b\\
    A=I_n, b=&[f(t_1), \ldots, f(t_n)]^T
\end{align*}
得到$x=[a_1,\ldots,a_n]^T$，故插值函数$P(x)=\sum_{i = 1}^{n} a_iB_i^1(x)$。

当order=3时，即三次样条插值。需要额外插入$t_{-2}\ldots t_{0}$和$t_{n+1}\ldots t_{n+3}$，根据
等距添加可以得到额外点。并且由于边界条件不同，A和b也不同：\\
$\bullet$way=1
\begin{equation*}
    A=\left(\begin{array}{ccccc}
        a_{1,1} & -a_{1,1}-a_{1,3} & a_{1,3} & &\\
        B_{-1}^3(t_1) & B_0^3(t_1) & B_1^3(t_1) & &\\
        & \ddots & \ddots & \ddots &\\
        & & B_{n-2}^3(t_n) & B_{n-1}^3(t_n) & B_{n}^3(t_n)\\
        & & a_{n+2,n} & -a_{n+2,n}-a_{n+2,n+2} & a_{n+2,n+2}\\
    \end{array}\right)
\end{equation*}
其中
\begin{equation*}
    a_{1,1}=-\frac{3B_0^2(t_1)}{t_2-t_{-1}}\qquad
    a_{1,3}=\frac{3B_1^2(t_1)}{t_3-t_0}
\end{equation*}
\begin{equation*}
    a_{n+2,n}=-\frac{3B_{n-1}^2(t_n)}{t_{n+1}-t_{n-2}}\qquad
    a_{n+2,n+2}=\frac{3B_n^2(t_n)}{t_{n+2}-t_{n-1}}
\end{equation*}
\begin{equation*}
    b=\left(\begin{array}{c}
        f'(t_1)\\
        f(t_1)\\
        \vdots \\
        f(t_n)\\
        f'(t_n)\\
    \end{array}\right)
\end{equation*}\\
$\bullet$way=2
\begin{equation*}
    A=\left(\begin{array}{ccccc}
        a_{1,1} & -a_{1,1}-a_{1,3} & a_{1,3} & &\\
        B_{-1}^3(t_1) & B_0^3(t_1) & B_1^3(t_1) & &\\
        & \ddots & \ddots & \ddots &\\
        & & B_{n-2}^3(t_n) & B_{n-1}^3(t_n) & B_{n}^3(t_n)\\
        & & a_{n+2,n} & -a_{n+2,n}-a_{n+2,n+2} & a_{n+2,n+2}\\
    \end{array}\right)
\end{equation*}
其中
\begin{equation*}
    a_{1,1}=\frac{6}{(t_2-t_{-1})(t_2-t_0)}\qquad
    a_{1,3}=\frac{6}{(t_3-t_{0})(t_2-t_0)}
\end{equation*}
\begin{equation*}
    a_{n+2,n}=\frac{6}{(t_{n+1}-t_{n-2})(t_{n+1}-t_{n-1})}\qquad
    a_{n+2,n+2}=\frac{6}{(t_{n+2}-t_{n-1})(t_{n+1}-t_{n-1})}
\end{equation*}
\begin{equation*}
    b=\left(\begin{array}{c}
        f''(t_1)\\
        f(t_1)\\
        \vdots \\
        f(t_n)\\
        f''(t_n)\\
    \end{array}\right)
\end{equation*}\\
$\bullet$way=3\\
此时的A与way=2时A一样。
\begin{equation*}
    b=\left(\begin{array}{c}
        0\\
        f(t_1)\\
        \vdots \\
        f(t_n)\\
        0\\
    \end{array}\right)
\end{equation*}\\
因此根据插值点，可以解线性方程组$Ax=b$，从而得到唯一的插值函数，进一步利用operator可以返回
插值函数在任一点的值。
\section{课后题实验结果分析}
\subsection{A}
在A.cpp文件中，对函数f分别进行了阶为3的ppForm插值和B\_spine插值，并且每种插值有针对三种不同边界条件
进行测试。并在区间内取100个等距点对插值函数进行测试，将插值函数的值输出到A\_opint.txt文件中。将每个
N的子区间中点处的插值误差向量的最大范数输出到A\_error.txt文件中。在A.py中将每种插值结果以及误差
可视化展示，生成图像如下：\\
$\bullet$ppForm:\\
\begin{figure}[H]
    \centering
    \subfloat[way=1]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_ppForm1.png}
    }
    \subfloat[way=2]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_ppForm2.png}
    }
\end{figure}
\begin{figure}[H]
    \centering
    \subfloat[way=3]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_ppForm3.png}
    }
    \subfloat[way=2]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_ppForm_error.png}
    }
\end{figure}
$\bullet$B\_spline\\
\begin{figure}[H]
    \centering
    \subfloat[way=1]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_B_spline1.png}
    }
    \subfloat[way=2]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_B_spline2.png}
    }
\end{figure}
\begin{figure}[H]
    \centering
    \subfloat[way=3]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_B_spline3.png}
    }
    \subfloat[way=2]
    {
        \includegraphics[width=0.5\textwidth]{figure//A_B_spline_error.png}
    }
\end{figure}
由图像可以看出，随着插值点的增加，插值函数与原函数f之间的图像越来越重合，插值函数值与原函数值
之间的误差越来越小，收敛度越来越高。并且无论哪种插值哪种边界条件，都具有很好的收敛效果，并且在
两端并未产生Runge现象。

下面针对收敛阶进行分析，由于题目所给N太少，导致无法明显看出收敛阶，因此增加三个N，使每次
插值点小区间个数分别为前一次的二倍。设$C_i$为第i次N插值的最大误差，则收敛阶p为
\begin{equation*}
    p=\lim_{i \to \infty}  \frac{C_i}{C_{i+1}}
\end{equation*}
这里为了方便观察，使比值取自然对数，以ppForm插值，way=1的情况为例，可以得到下图：
\begin{figure}[H]
    \centering
    \subfloat[degree of convergence]
    {
        \includegraphics[scale=0.45]{figure//A_convergence.png}
    }
\end{figure}
由图像可知，随着插值点的增加，误差比趋于2.8（以e为底），故收敛阶为2.8。
\subsection{B}
问题B在头文件spline.h中实现，class B\_spline可以分别求出不同阶数的B样条基函数，再通过class
B\_spline\_interpolation可以实现不同边界条件的三阶B样条插值，并且由于线性方程组的唯一解可以
得出插值函数的唯一性。
\subsection{C}
针对Theorem 3.57，可以采用way=1、order=3的B样条插值，对于Corollay 3.58，采用order=2的B样条插值。
同时分别取1000和900个值代入插值函数，并将插值函数值分别输出到.txt文件，并利用python可画出图像。但是
由于Corollary 3.58所采用B样条基函数为以整数为区间端点，而插值点并非整数，因此我尝试再单独写一个特殊
B样条整数基函数，无奈能力所限，写了一天最终还是没有成功，因此图像展示的二次样条图像与原函数略有不同。
为了可以继续做下去，我增加了对一次线性样条的分析，故图像如下：
\begin{figure}[H]
    \centering
    \subfloat[the polynomials against the exact function]
    {
        \includegraphics[scale=0.45]{figure//C.png}
    }
\end{figure}
\subsection{D}
根据插值误差公式$E_S(x)=|S(x)-f(x)|$,在题目所给点对C题中的三次样条和一次样条分别计算插值误差，
并将误差结果存入D\_error.txt文件中，其输出结果如下图：
\begin{figure}[H]
    \centering
    \subfloat[插值误差结果]
    {
        \includegraphics[scale=0.45]{figure//D_error1.png}
    }
\end{figure}
可以看出有部分x处的误差接近机器精度，即为0，并且这些x为插值节点。主要原因是样条插值要求
插值节点处插值函数值与原函数值相等，故此处误差接近机器精度。为了判断两个不同基数插值的精确度差异，
将插值误差函数可视化，如下图：
\begin{figure}[H]
    \centering
    \subfloat[$|S(x)-f(x)|$]
    {
        \includegraphics[scale=0.45]{figure//D_error.png}
    }
\end{figure}
可以明显看出，三次样条插值误差小于一次样条插值误差，显然三次样条插值更加精确。
\subsection{E}
根据题目所给心形曲线方程可以得到：
\begin{align*}
    x^2+(\frac{3}{2}y&-\sqrt{|x|})^2 = 3\\
    y = \pm \frac{2}{3}\sqrt{3-x^2}&+\frac{2}{3}\sqrt{|x|} \qquad x^2 \leq 3 
\end{align*}

因为心形曲线为二维，故需要两次样条插值。由于该曲线较为复杂，很难给出待插值函数的解析表达式，导数
也不是那么容易计算，而自然样条插值的边界条件不涉及导数，故作为心形曲线插值的边界条件最为合适。另外，
不管是由于该图像为曲线还是前面题目结果所得，三次样条插值比一次样条插值精确度显然更高。综合来看，
采用三次样条插值，以自然样条作为边界条件最为合适。

由于需要两次样条插值，并且需要利用累积弦长进行计算。故在E1.cpp中写了一个计算离散点值的class类。
最终对心形函数进行区间划分分别作出图像，可以得到如下图：
\begin{figure}[H]
    \centering
    \subfloat[n=10]
    {
        \includegraphics[width=0.5\textwidth]{figure//E1.png}
    }
    \subfloat[n=40]
    {
        \includegraphics[width=0.5\textwidth]{figure//E2.png}
    }
\end{figure}
\begin{figure}[H]
    \centering
    \subfloat[n=160]
    {
        \includegraphics[width=0.5\textwidth]{figure//E3.png}
    }
    \subfloat[对比图]
    {
        \includegraphics[width=0.5\textwidth]{figure//E.png}
    }
\end{figure}
显然可以看出，随着n的增大，心形曲线更加圆滑。
\end{document}